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ABSTRACT 

Cool cores of galaxy clusters are thought to be heated by low-power active galac- 
tic nuclei (AGN), whose accretion is regulated by feedback. However, the interaction 
between the hot gas ejected by the AGN and the ambient intracluster medium is ex- 
tremely difficult to simulate as it involves a wide range of spatial scales and gas that is 
Rayleigh- Taylor (RT) unstable. Here we present a series of three-dimensional hydrody- 
namical simulations of a self-regulating AGN in a galaxy cluster. Our adaptive-mesh 
simulations include prescriptions for radiative cooling, AGN heating and a subgrid 
model for RT-driven turbulence, which is crucial to simulate this evolution. AGN 
heating is taken to be proportional to the rest-mass energy that is accreted onto the 
central region of the cluster. For a wide range of feedback efficiencies, the cluster reg- 
ulates itself for at least several 10 9 years. Heating balances cooling through a string of 
outbursts with typical recurrence times of around 80 Myrs, a timescale that depends 
only on global cluster properties. Under certain conditions wc find central dips in the 
metallicity of the intracluster medium. Provided the sub-grid model used here captures 
all its key properties, turbulence plays an essential role in the AGN self-regulation in 
cluster cores. 
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1 INTRODUCTION 



High-resolution observations from the Chandra and XMM- 
Newton X-ray telescopes have revolutionised our under- 
standing of the hot-intracluster medium (ICM). While a 
large number of galaxy clusters have strong peaks in their 
central X-ray surface brightness distributions, indicating 
that their central gas is cooling rapidly, detailed spectra of 
such cool-core clu sters show that this gas fails to cool be- 
low fts 1 keV (e.g. IPeterson et allkoOll ; IRaffertv et"ai1l2006l : 
iMcNamara fe Nulsertl2007h . This deficit means that radia- 
tive cooling must be balanced by an unknown energy source. 

Currently, the most successful model for achieving 
this balance relies on heating by ou tflows from an AGN , 

|2000| ; 



hosted by the cent r al cD galaxy (e.g. Churazo v et al 



iBriiggen fc Kaiseil 120021 : iMagliocchetti fc Briiggenl l2007h 
Galaxies within 0.2r2oo of the centre of groups and clusters 
show a boosted likeliho od of hosting large radio-lou d jets 
of AGN-driven material (|Burnslll990l : IBest et al.ll2005h . The 
energies from such outbursts are co mparable to those neede d 
to stop the gas from cooling {e.g. ISimionescu et af1 l2008). 
and the mean powe r of the outbursts is well-correlated with 
the radiated power (Birz an et al.| [2004). Furthermore, AGN 
feedback from the central cD increases in proportion to the 



cooling lu minosity, as expected in an operational feedback 
loop (e.g. iBirzan et al]|2004 IRaffertv et al.lbOOfj '). 

Yet despite these many clues, the nature of the feed- 
back mechanism is still not understood. In clusters, energy 
coming from a region less than a parsec across has to regu- 
late a volume of gas that is about a few hundred kpc across. 
This regulation must maintain the balance between heat- 
ing and coo ling for a long tim e, at least since a redshift 
of z m 0.4 ijBauer et al.l 120051 ). Moreover, it must hold in 
systems wit h luminosities that sp an about three orders of 
magnitude (jChurazov et al.| [2005). Finally, the output by 
the AGN must be regulated by the physical conditions in 
the vicinity of the AGN, rather tha n the properties of the 
host galaxy in which it is contained (|Best et al.l 120071 ). 

In general, gas can make its way onto an AGN disc 
directly from the hot gas via Bondi-Hoyle accretio n (e.g. 
Allen et al.l 12006m or in the form of cold blobs (e.g. ISokeil 



2006). It is likely that the AGN feedback in clusters is 



caused by accretion of hot gas onto pr e-formed massive 
black holes in old ellip tical galaxies (e.g. iBest et al 1 120071 : 
lHardcastle et al.ll2007n . which produce low-power radio jets 
but relatively little emission at optical, ultraviolet and in- 
frared wavelengths. On the other hand, the fuelling by cold 
gas m ay be responsible for A GN with strong emission lines 
(e.g. Hardca stle et al.l [2007') but the supply of powerful 
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high-excitation sources does not have a direct connection 
to the hot phase. Most efforts to explain feedback have con- 
cent rated on the hot - accre tion scenario. 

[Ch urazov et all (|2002h considered a toy model of ICM 
feedback regulated by hot gas accretion, in which black 
hole growth was set by the classic Bondi formula, with 
M oc n e T~ 3 / 2 . For a ratio of specific heats of 5/3, this 
rate only depends on the entropy and is directly affected 
by heating and cooling. The ICM in the model responded 
to the heating by expanding, which lowered the gas density 
and accretion rate. As the gas radiated away energy, the at- 
mosphere contracted and the accretion rate rose again. Ap- 
plication of the results to the Virgo cluster suggested that 
accretion proceeded at approximately the Bondi rate down 
to a few gravitational radii and that most of the power was 
carried away by an outflow. Comparison of the accretion 
rate predicted by the Bondi formula with the luminosity of 
the cooling flow region in M87 suggested a high efficiency 
for transforming the rest mass of accreted material into me- 
chanical power of an outflow at the level of a few per cent. 
Recently, I Allen et"al1 (|2006l ) found evidence that the power 
of outbursts in giant elliptical galaxies scales with the Bondi 
accretion rate. 

However, the full simulation of this feedback cycle re- 
mains elusive. Direct simulation of a galaxy cluster from 
Mpc scales down to the parsec scale at which mass ac- 
cretes onto the central supermassive black hole is still 
out of reach even for the latest generation of adaptive- 
mesh refinement (AMR) codes. Hence, in most hydro- 
dynamical simulations of AGN feedback, the energy in- 
put into the ICM is set by hand, rather than computed 
from the conditions surrounding the central galaxy. Onh 



IVernaleo fc Reynolds! f2 006). Br ighenti fc Mathews! ([200' 



and ICattaneo fc Tevssierl (|2007h attempted hydrodynamic 
simulations of self-regulated accretion, and none of these 
studies was able to reproduce all the observed properties of 

cool -core clusters. 

IVernaleo fc Reynolds! l|2006l ) carried out a set of three- 
dimensional cluster simulations in which they took the ki- 
netic energy of the AGN jet to be proportional to the in- 
stantaneous mass accretion rate across an inner spherical 
boundary. They found that their simulations were incapable 
of balancing heating and cooling as the jet formed a low- 
density channel, which stopped the hot gas from coupling 
effectively to the surrounding ICM. A delay between the 
mass accretion rate and the response of the jet was intro- 
duced in further si mulations, but cou ld not alleviate this 
problem. However, iHeinz et alj l|2006l ) have shown that in 
more realistic initial conditions the ambient bulk motions 
of the ICM prevent the formation of such channels, though 
thei r simulations were not s elf-reg ulating. 

iBrighenti fc Mathews! |2006) simulated self-regulated 
mass-loaded jets that resemble disk winds. They presented 
two-dimensional hydrodynamical models where outflows 
were generated by assigning a fixed velocity to gas that 
flowed into a biconical source region at the cluster centre. 
The acceleration of gas in the source cones was triggered 
when gas flowed into the innermost zones. As the jet pro- 
ceeded, it entrained additional ambient gas and its mass flux 
increased. After many gigayears, the time-averaged gas tem- 
perature profile resembled observations and very little gas 
cooled below the virial temperature. 



However, feedback is at least partly mediated by low- 
power Fanaroff-Riley t ype I (FR I) sources a n d this has not 
been simulated before. ICattaneo fc Tevssierl (|2007h carried 
out three-dimensional simulations in which AGN jets were 
regulated via a Bondi-accretion rate measured at the cluster 
centre starting from hydrostatic initial conditions. For their 
setup, a very violent AGN phase (L ~ 10 45 erg/s) occurred 
lasting for about 250 million years, followed by a period of 
at least 4 gigayears during which the accretion rate of the 
black hole and the mass of the cold gas in the cluster core 
stayed constant with a luminosity of a few times 10 43 erg 
s _1 . This scenario does not explain the recurrent, low-power 
FR I sources that are observed in local cool-core clusters. 

The biggest obstacle to these simulations is that it is 
still unclear how the mechanical power of the AGN is trans- 
formed into thermal energy and distributed throughout the 
cluster. iNulsen et al.l l|2007h pointed out that theoretically 
one would expect the energy to be dissipated locally. As the 
bubble rises, they argue, the ICM falls inward around it to 
fill the space it occupied, converting gravitational potential 
energy to kinetic energy. Details of the flow depend on the 
viscosity, which is poorly determined. If it is high, the flow 
is laminar and the kinetic energy is dissipated as heat over 
a region comparable in size to the cavity. If the viscosity is 
low, the Reynolds number is high and the flow would form 
a turbulent region of a similar size to the cavity. Turbu- 
lent kinetic energy is dissipated in the turnover time of the 
largest eddies, so that the dissipation time fdiss ~ '■fc/tWb, 
where rb is the radius of the bubble and Uturb is the speed 
of the eddy, which is comparable to the speed of the cav- 
ity. Since tdiss^b ~ f&, much of the turbulent kinetic energy 
is dissipated in a region of comparable size to the cavity. 
Thus, regardless of the viscosity, the kinetic energy created 
by cavity motion is deposited locally. Other authors have 
pointed out the imp ortance of sound waves and weak shocks 
in dissipating heat dFabian et"al"1l2003l ; iBriiggen et al.ll2007l ; 
ISanders fc FabianlbOOTT 2008), and considered the possibil- 
ity that AGN bubbles are stabilized by magnetic fields (De 
Young 2003; Bri iggen fc Kaiser 2001; Rus zkowski et al. 2007) 
or massive iets ^Sternberg fc SokeJl2008l ). 

In lScannapieco fc Briiggenl l|2008h (SB08 hereafter), we 
showed that even in the absence of magnetic effects, the evo- 
lution of AGN bubbles cannot presently be captured reliably 
by direct simulations. While pure hydrodynamic simulations 
indicate that AGN bubbles are disrupted into resolution- 
dependent pockets of underdense gas, proper modeling of 
subgrid turbulence indicates that this a poor approxima- 
tion to a turbulent cascade that continues far beyond cur- 
rent resolution li mits. Adding a subgrid model of turbu- 
lence and mixing (|Dimonte fc Tiptonll20od ) to the adaptive- 
mesh hydrodynamic code, FLASH3, we found that Rayleigh- 
Taylor (RT) instabilities act to effectively mix the heated re- 
gion with its surroundings, while at the same time preserv- 
ing it as a coherent structu re, consistent with observations. 
ICattaneo fc Tevssierl (|2007h also pointed out the importance 
of turbulence in mediating the jet power to the accretion 
mechanism, although they were unable to resolve it. 

In this paper, we explore how a subgrid model for tur- 
bulence provides the missing piece necessary to construct a 
full numerical simulation of a self-regulating jet at the centre 
of cool-core clusters. Furthermore, with no fine tuning, such 
simulations naturally explain both the mechanical power of 
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AGN, and their duty cycles. The structure of the paper is 
as follows: In §2, we describe our method including the al- 
gorithm and the setup, and in §3 we present our results. A 
discussion is given in §4. 



2 METHOD 

2.1 Simulation and Turbulence Modeling 

Out simulations we re performed with FLASH version 3.0 
jFrvxell et al.ll2000t ). a multidimensional adaptive mesh re- 
finement hydrodynamics code. It solves the Riemann prob- 
lem on a Cartesian grid using a directi onally-split Piecewise- 
Para bolic Method (PPM) solver (|Colella fc Woodward] 
1 19841 ) . While the direct simulation of turbulence is extremely 
challenging, com putationally expen sive, and dependent on 
resolution (e.g. iGlimm et alj 120011 ) . its behaviour can be 
approximated to a good de gree of accuracy by adopt ing a 
sub-grid approach. Recently. iDimonte fc Tiptonl (120061 ). de- 
scribed a sub-grid model that is especially suited to captur- 
ing the buoyancy-driven, turbulent evolution of AGN bub- 
bles. The model captures the self-similar growth of the RT 
and Richtmyer-Meshkov (RM) instabilities by augmenting 
the mean hydrodynamics equations with evolution equa- 
tions for the turbulent kinetic energy per unit mass and the 
scale length of the dominant eddies. The equations are based 
on buoyancy-drag models for RT and RM flows, but con- 
structed with local parameters so that they can be applied to 
multi-dimensional flows with multiple materials. The model 
is self-similar, conserves energy, preserves Galilean invari- 
ance, and works in the presence of shocks. 

The mean flow fluid equations in this case are given by 



this as 



dp dpuj 
dt dxi 



0. 



dpiii dpuiiij 
dt dxj 



dP 9$ _ pK 

+ a h P~x~ = Cp "5 — ' 

oxi oxi axi 



(1) 



(2) 



dpE dpEuj dP 



dt 



+ 



dxj 



+ 



dxi 



JL 

dx 



V N E dxj J 



where t and x are time and position variables, p(x, t) is the 
average density field, fij(x, t) is the mass-averaged mean- 
flow velocity field in the i direction, P(x, i) is the mean 
pressure, $(x,t) is the gravitational potential, and E(x.,t) 
is the mean internal energy per unit mass. The terms on 
the right hand side of these equations couple the subgrid 
quantities to the mean flow. These three terms model: (i) the 
excess "pressure" arising from the turbulent kinetic energy 
per unit mass K, which is scaled by a constant Cp; (ii) 
turbulent mixing, which is modeled as a turbulent viscosity 
pt, scaled by a constant Ne; and (iii) an explicit source 
term, Sk that contains both RT and RM contributions, and 
removes energy from the mean flow and moves it into the 
turbulent field. This is given by 



S K = pV 



CE-AiQi ~ Cd 



V 2 



(3) 



where V = V2K is the average turbulent velocity, L is the 
turbulent eddy scale, gi = — (1/ p)dP/dxi is the gravita- 
tional acceleration, the coefficients Cb and Cd are fit to 
experiments (Cb = 0.84 and Cd = 1.25) and Ai is the At- 
wood number in the i direction. In the code, we determine 



At = 



P++P- 



dp 



p + L\dp/dxi\ dxi ' 



(4) 



where Ca is a constant and p_ and p+ are the densities 
on the back and front boundaries of the cell in the i direc- 
tion. This source term causes turbulence to decay away at a 
characteristic time scale of w L/V in the absence of external 
driving, while the growth term causes turbulent velocity to 
grow at a rate V ~ g, whenever the gravitational accelera- 
tion opposes the density gradient. 

The turbulence quantities that appear in these equa- 
tions are calculated from evolution equations for L and K. 
The eddy scale L must be evolved because the buoyancy- 
driven RT and RM instabilities depend on the eddy size, 
which is expected to grow self-similarly. Simple equations 
that include diffusion, production, and compression arc 
given by 



dpL dpLHj 

dt dxj 

and 

dpK dpKuj 



d_fjH_dl/. 
dxi V N L dx, 1 ' 



CapL^- , (5) 



dt 



dxj 



d (jH_dK\ 

dxj \ Nk dxj J 



CppK^ + Sk, (6) 



where Nk, Nl, Co, and Cp are con stants (Nk = 1-0, 
N L = 0.5, C c = 1/3 and C P = 2/3). See lDimonte fc Tiptonl 
(2006) for details on how these constants are determined. In 
the L equation the three terms on the right hand side of the 
equation represent, respectively: turbulent diffusion, growth 
of eddies through turbulent motions, and growth of eddies 
due to motions in the mean flow. In the K equation the three 
terms on the right hand side represent, respectively: turbu- 
lent diffusion, the work associated with the turbulent stress, 
and the source term Sk, which also appears in the internal 
energy equation to conserve energy. Finally, the turbulent 
viscosity is calculated as 



Pt = Cy,pL\ 



(7) 



where C M is a constant. Further details of this model and 
our implementation are given in SB08. 

All our simulations are performed in a cubic three- 
dimensional region with reflecting boundaries that is 680 
kpc on a side, although the cool-core region in which we are 
most interested lies within < 100 kpc of the center. For our 
grid, we chose a block size of 8 3 zones and an unrefined root 
grid with 8 3 blocks, for a native resolution of 10.6 kpc. The 
refinement criteria are the standard density and pressure cri- 
teria, and we allow for 4 levels of refinement beyond the base 
grid, corresponding to a minimum cell size of 0.66 kpc, and 
an effective grid of 1024 3 zones. In all zones K is initialized 
to 10~ 10 (p/p)/2. and L is initialized to 10~ 5 y^p/py/l/Gp. 

2.2 Cluster Profile 

For our o verall cluster p r ofile, we adopted the model de- 
scribed in iRoedieer et"atl \200H ). which was constructed to 
reproduce the properties of the brightest X-ray cluster A426 
(Perseus) that has been studied extensively with Chan- 
dra and XMM-Newton. In this case, the electron density 
n c and temperature T c p rofiles are based on the depro- 
jected XMM-Newton data (jChurazov et alj |2003) which are 



© 2009 RAS, MNRAS 000.1111171 



4 M. Briiggen & E. Scannapieco 



in broad agreement with t he Chandra data [1 Schmidt et al.l 
|2002| ; ISanders et aJj|2003 ); namely: 



4.6 x 10" 



+ 



4.8 x 10" 



[l+(£) 2 ] 18 [1 + (w) 2 ] ai 



and 



T c = 7x 



[i + (tt) 3 ] 

[2-3 + (£) 3 ] 



keV, 



(8) 



(9) 



where r is measured in kpc. Furthermore, the hydrogen 
number density was assumed to be related to the electron 
number density as tig = n c /1.2. The static, spherically- 
symmetric gravitational potential, $(-??), was set such that 
the cluster was in hydrostatic equilibrium. 

As in Scannapieco & Briiggen (2008), we computed the 
radiative losses in each cell in the optically-thin limit, with 
an emissivity e = A(T)n_Hn e where hh and n e are the 
number density of hydrogen and electrons, respectively, and 
the cooling function A(T) descri bes radi a tive l osses from 
the optically thin plasma, as in ISarazml (|l986l ) (see also 
iRavmond et al.lll976l ; iPeres et allll982ft . 

The metal injection rate in the central galaxy was as- 
sumed to be proportional to the light distribution, and mod- 
eled with a Hernquist profile given by 



Pmetal(r, t) 



A/mcta lO a 



(r + a) 



(10) 



with a scale radius of a = 10 kpc and a total metal injection 
rate of M me taio = 9-5 x 10 4 Af© Myr -1 . This is the integral 
mass injection rate for all elements heavier than helium. The 
sources of metals ar e supernovae la and stellar mass loss 
jRebusco et alj|2005l ). 

In order to be able to trace the metal distribution, we 
utilize a mass scalar to represent the local metal fraction in 
each cell, F = p mG t a i/p. Hence, the quantity Fp gives the 
local metal density, p me tai, which has a continuity equation 
including the metal source, given by 



dpF dpFiij _ _9_ f JM_9F^ + p meta i 
dt dxj 



dij \ Nf dxj J p 



(11) 



where Nf is a dimensionless constant. Furthermore, we as- 
sumed that the metal fraction is small at all times. Hence, 
we could neglect p mo tai as a source term in the continuity 
equation for the gas density. Also, the mass fraction can be 
scaled to any value in post-processing. 

2.3 AGN Feedback 

Bubbles in the ICM are thought to be inflated by a pair 
of ambipolar jets from the central AGN that inject energy 
into small regions at their terminal points, which expand 
until the y reach pressure equilib rium with the surrounding 
medium piandford fc Reeslll974r i. The result is a pair of un- 
derdense, hot bubbles on opposite sides of the cluster centre. 
In order to produce such bubbles, we injected energy into 
two small spheres of radius rabble, each a distance rbubbie 
from the centre. 

Unlike in our previous simulations, the energy input 
rate was not predetermined, but instead calculated from the 
instantaneous conditions near the centre of the cluster. In 
particular, we considered a model in which a fixed fraction of 
the gas within the central 3 kpc of the cluster accretes onto 



Table 1. Run Parameters 



Run 


Subgrid 


-^bubble 


r bubble 


Name 


Model 


(ergs) 


(kpc) 


N5-10 


No 


5 x 10 59 


10 


D5-10 


Yes 


5 x 10 59 


10 


D5-8 


Yes 


5 x 10 59 


8 


D5-12 


Yes 


5 x 10 59 


12 


D2-10 


Yes 


2 x 10 59 


10 


D10-10 


Yes 


1 x 10 60 


10 


D20-10 


Yes 


2 x 10 60 


10 



the central supermassive black hole within a cooling time, 
and a fixed fraction of the me 2 rest mass energy of this ac- 
creted gas is returned to the ICM by increasing the pressure 
within t he bubble regions. Th i s is si milar to the approach 
taken bv lVernaleo fc Reynold si (|2006n . who calculated AGN 
feedback based on the mass inflow rate along the 10 kpc 
inner boundary of a spherical grid. Here our choice of the 
volume within 3 kpc of the cluster centre is taken to sample 
the conditions as near to the central supermassive black hole 
as possible, while still averaging over enough cells to avoid 
numerical effects. 

The efficiency of gas accretion in powering the bubbles 
was chosen by scaling to typical numbers observed in galaxy 
clusters. We took the rate of change of the pressure within 
the bubble regions to be given by 



-Pi 



, , 47T 3 
nj W - ^ -r bubblc 



E b 



75Myrs Af gas ,3kpc(t) 

tcool,3kpc(t) Af gaSi 3kpc(0) 



(12) 



where t C ool,3kpc(i) is calculated by dividing the total inter- 
nal energy of the gas within the central 3 kpc by the total 
radiation rate in this regions, and M ga ,s,3k P c (i) is the to- 
tal gas within the central 3 kpc, and -Ebubbic ~ 5 x 10 59 
ergs corresponds to the typical energies necessary to in- 
flate the X-ray cav ities observed in cool-core clusters (e.g. 
iForman et~ai1l2007i ) . which are generated approximately ev- 
ery 75 Myrs. This prescription corresponds to an overall 
efficiency of 3.2 x 10" 3 £ b ubbi e /10 60 ergs of the mc 2 energy 
of the mass cooling rate within the central 3 kpc. Note that 
choosing a typical value of -Ebubbic ~ 5 x 10 59 ergs would 
mean that our energy input is only about 4 x 10" 5 of the 
rest mass energy of the mass cooling rate within the central 
10 kpc, and thus our models generally have energy input 
levels similar to that considered by Vernaleo and Reynolds 
(2006). 

Finally, we started all runs with an initial pair of bub- 
bles overpressured by a factor of 



-Pbubble(O)— r bubblc 



(13) 



so as to create some initial "seed" turbulence, to help in- 
creasing mixing at early times and speed the simulations 
towards their eventual quasi-steady states. The choices of 
^bubble and r b u bb i c for each of the runs carried out in this 
study are given in Table [T] 
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Figure 1. Bubble power and integrated cooling rates for D5-10 (top left), N5-10 (top right), D5-8 (bottom left) and D5-12 (bottom 
right). The solid blue lines shows the cooling losses in the entire simulation volume, the dashed blue line show the cooling losses in 
the central 100 kpc of the cluster and the red line shows the energy put into bubbles. The lines for our fiducial run without a subgrid 
turbulence model (named N5-10) stop at about 450 million year because this is where catastrophic cooling in the centre halted the 
simulation. 



3 RESULTS 

3.1 Self-regulating Oscillations 

As a control case, we first present the results of a model 
without subgrid turbulence (5N-10) in which we have chosen 
^bubble ~ 5 x 10 59 ergs and rbubbie = 10 kpc. The evolution 
of the cooling and AGN heating in this case is summarised in 
the top left panel in Fig. [T] Here we see that the bubbles do 
not couple to the central region that determines the activity 
of the AGN. Instead, cool gas accumulates in the centre, 
increasing the AGN output as shown in this figure, as well 



as leading to the cold dense central regions shown in Fig. [2] 
This cooling eventually leads to a drastic outburst of the 
AGN, which is apparent in a sharp upturn in the heating 
rate in Fig. [T] and the large empty cavities in the third row 
in Fig. [2] However this heating does not stop the flow of gas 
onto the AGN along the midplane. Instead cooling increases 
catastrophically and we are forced to stop the simulation. 

On the other hand, when the subgrid model for the tur- 
bulence is switched on, as in run 5D-10, the bubbles couple 
much more effectively to the central region. As seen in the 
second and third columns of Fig. [5] the cold gas that forms 
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20 40 60 80 100 120 140 20 40 60 80 100 120 140 20 40 60 80 100 120 140 

Mpc Mpc Mpc 




-26.0 -25.8 -25.5 -25.2 -25.0 -24.8 -24.5 0.6 0.7 0.8 0.9 1.1 1.2 1.3 0.9 1.3 1.7 2.1 2.6 3.0 3.4 



Figure 3. Slices of log(p/gcm 3 ) (1st column), log sound crossing time in Myrs (2nd column), and log turbulent diffusion time in Myrs 
(3rd column), as the fiducial subgrid run (D5-10) executes a full ~ 80 Myr cycle as described in the text. From top to bottom t = 600, 
620, 640, 660 and 680 Myrs (all from D5-10). 
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Figure 4. Volume rendering of the temperature in our fiducial run D5-10 at t= 630 Myrs. The blue ring shows the cool gas accreting 
onto the centre; the red and yellow blobs represent the hot gas ejected by the AGN. On the left side one can see an older bubble from 
an earlier outburst. 



in a ring around the ICM, shown in white, is slowly eroded 
and the AGN activity subsides. The turbulent diffusion or 
turn-over time, L/V, is plotted in the third column of Fig. [2] 
It shows that near the ring of cool gas in the centre of the 
cluster this time is of the order of 80 Myrs. This turbulence 
mixes in hot material near the centre of the cluster, heating 
the cool inflowing ring of intracluster gas and stopping the 
AGN outburst. 

Thus instead of cooling catastrophically, the cluster in 
the simulation with subgrid turbulence begins to execute a 
series of oscillations as indicated in Fig. [3] The top row of 
this plot shows the configuration at 600 Myrs, shortly after 
an outburst of AGN activity. From this time to 640 Myrs, 
turbulence decays away at the turn-over time scale, L/V, 
and mixing near the cluster center becomes progressively 
less efficient. This leads to an increased level of accretion 
near the core, as more and more cold gas makes its way 
onto the AGN. The result is a new burst of AGN activity, 
which drives a bubble on the order of a sound crossing time. 
The RT-unstable bubble leads to a rise of the turbulent vis- 
cosity and mixing, which quickly quenches accretion when 
the turbulent length scale grows to be of the order of the 
scale of the accretion flow. At this point the AGN remains 
relatively quiescent until the turbulence decays again on a 
time of order L/V , repeating the cycle as shown in the bot- 
tom panels of Fig. [3] This interaction of the AGN-heated 
regions with the cool inflowing gas is also illustrated in the 
volume-rendered image shown in Fig. [3] 

It is important to point out that we did not tweak 



the parameters of the subgrid model to achieve the self- 
regulating cycle. Instead, these parameters are determined 
by experimental studies of the RT instability and are th e 
same as those employed in IScannapieco fc~B riiggcn (2008). 
Note also that the time scale for this cycle is not the sound 
crossing time l/c s for the central region, which is of the order 
of 10 Myrs. Instead the period between the episodes of AGN 
fueling is determined by the time it takes for the turbulence 
to decay away after it has begun to mix the cold gas in the 
cluster centre. This is given by 



td 



uty 



I/Vtu 



rb ■ 



(14) 



where I is the distance between the accretion region of the 
AGN and the heating region. In our fiducial simulation, 
I = ^bubble = 10 kpc and is clearly resolved, iwb is a typical 
turbulent velocity, which we assume grows in a dynamical 
time given by «tmb ~ gt ~ gl/c s , and c s is the sound speed. 
Because the cluster is in hydrostatic balance, the gravita- 
tional acceleration, g can be written as 



2 1 dp _ 2 i 
9 = c sz~ =Cs/ro- 



' p dr 

This leads to the turbulent velocity 



and 



"turb = gi/Cs ~ C s , 

ro 



/■() 



tduty ^ 



(15) 



(16) 



(17) 
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Figure 5. Bubble power and integrated cooling rates for D5-10 (top left), D2-10 (top right), D10-10 (bottom left) and D20-10 (bottom 
right). Meaning of lines as in Fig. 1. 



According to this estimate the turbulent speed is given by 
c 3 l/r , which is roughly 700 km s" 1 (10 kpc/ 60 kpc) ~ 100 
km s ~ , which is indeed what we find in our simulations. 
This means that the size of the bubbles does not enter the 
expression for the duty cycle as would be expected in a cycle 
regulated by laminar flow. It is the properties of the cluster 
itself, rather than the jet physics of the central radio source 
that are setting the recurrence time of the bubbles. Thus 
the size of the bubbles themselves should not affect the re- 
currence time of the AGN that we see in Fig. [T] - [3j For the 
cluster simulated here, the central sound speed is around 
700 km s~\ the scale height of the cluster, rrj, is around 60 
kpc so the duty cycle is 60 kpc / 700 km s" 1 w 80 Myrs. 

To test this hypothesis we have varied the geometry 



of the injection region, leading to the heating and cooling 
evolution shown in the lower panels of Fig. [T] As in run 
D5-10, the time between two subsequent outbursts for runs 
D5-8 and D5-12 is roughly 80 Myrs. Also, in all cases the 
activity goes down when the cool core has been destroyed 
by recurrent AGN activity. Again, this large overall time 
scale, which is w 1 Gyr, is not dependent on the particulars 
of bubble injection. Rather it is the time it takes for the 
turbulence to dissipate through the cool core of the clusters, 
which is given by 



■Rcc/wturb ~ 100kpc/100km/s w 1 Gyr 



(18) 



After this the core has been heated sufficiently to stop the 
gas inflow. It then takes a central cooling time of about 1 
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Figure 6. Slices of the logarithm of the density at 1 Gyr. Top from left to right: D5-10, D5-8, D5-12 Bottom from left to right: D2-10, 
D10-10, D20-10. 



Gyr for the cool core to reform and the intermittent AGN 
activity to resume. 

In Fig. we vary the efficiency with which inflowing 
rest mass energy is converted into thermal energy deposited 
by the AGN. The heating/cooling curves for runs D2-10, D5- 
10 and D10-10 look remarkably similar. Even with different 
efficiencies, the AGN easily regulates itself and the resulting 
duty cycles are largely unchanged. 

In fact, the only run that shows a noticeably different 
evolution history is the high energy D20-10 case. Although 
this case shows a strong AGN event after about 400 Mrs as 
in the other runs, unlike the other runs the energy injection 
remains relatively weak at late times and varies irregularly. 
So, provided the energy is injected sufficiently far from the 
region that determines the fuelling of the AGN and provided 
the heating rate is not so high as to destroy the cool core im- 
mediately, the duty cycle is independent of the exact choice 
of parameters. 



3.2 Overall Cluster Properties 

Slices of the density through the cluster centre at t = 1 Gyr 
and t = 1.5 Gyr are shown in Fig. [6] -[7] all plotted on the 
same colour scale. In these figures, the top rows show results 
from runs with varying bubble sizes at a fixed feedback ef- 
ficiency, and the bottom rows show results from runs with 
varying efficiencies at a fixed bubble size. The t = 1 Gyr 
slices shown in Fig. [S] correspond to a time shortly after a 



major outburst has taken place in each of the runs. One can 
see large plumes streaming from the centre in all cases, al- 
though these are somewhat weaker in the D5-8 run than in 
the other runs. In general, the size of the energy injection re- 
gion has a stronger impact on the appearance of the cluster 
than the choice of feedback efficiency, as one would expect 
for a self-regulating cycle, in which feedback increases only 
to the level necessary to balance cooling. In fact, for efficien- 
cies that vary by a factor of 10, the density distributions look 
remarkably similar. In the t — 1.5 Gyr slices, shown in Fig. 
[7] one only sees traces of the AGN activity in the form of 
bubbles of radius m 10 kpc on either side of the centre. Again 
these bubbles are clearly seen in all runs, and the variance 
between runs is more strongly dependent on geometry than 
feedback efficiency. 

X-ray maps corresponding to the same times as in the 
density slices are shown in Figs. [8] and [9] For simplicity we 
construct these by projecting the full cooling emissivity as 
computed from e = n_r/n e A(T) as described in §2.2. Because 
the emissivity is proportional to the square of the density, 
the X-ray images are dominated by the emission from the 
cluster centre, and demonstrate the same overall pe aked dis- 
tribu tion observed in nearby cool-core clusters (e.g. ISarazinl 
1986). On closer inspection, one can also see substructure 
and signatures of the AGN-blown cavities. In the 1 Gyr pro- 
files, for example, many runs show clear evidence of circular 
cavities that are surrounded by brighter emission, especially 
on the sides. This is particularly clear in run D5-12, which at 
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Figure 7. Slices of the logarithm of the density at 1.5 Gyr. Top from left to right: D5-10, D5-8, D5-12. Bottom from left to right: D2-10, 
D10-10, D20-10. 



1 Gyr contains a large X-shaped structure that resembles the 
X-ray image of the Virgo elliptical M84 (jFinoguenov et all 
2008). In all the runs, and unlike in all previous attempts 
to simulate self-regulating AGN, the appearance of our self- 
regulated AGN in X-rays resembles observations of nearby 
clusters (e.g. Fabian et al. 20031; McNamara et al. 20051 ; 



iNulsen et al.ll2005l : iForman et al.ll2007l : iKraft et alJl2007t >. 

The radial profiles of average density, emission- weighted 
temperature, and entropy are shown in Figs. 1101 - 1111 Here 
the density is given by an average over the solid angle, fl, as 

1 

47T . 



p(r) 



dO, p(r, fl), 



(19) 



the emission-weighted temperature is given by 

""*>= fd f%? < 2o » 

J ail e(r, 11) 

and the radial entropy profiles is calculated from the 
emission- weighted temperature according to 

S(r) = k B T EW (r)n c (r)- 2/s . 

In these figures, the right columns show the 
varying efficiencies, D2-10, D5-10 and D10-10, and the left 
columns show the runs with different injection geometries, 
D5-8, D5-10 and D5-12. 

Across the runs and at both time steps, the density con- 
tours deviate very little from the initial profile, as shown by 
the thick solid lines. On the other hand, there are notice- 
able changes in the temperature, which are also reflected 



(21) 
runs with 



in the entropy profiles. These differences are largest within 
50 kpc, and at these distances, the profiles at 1 Gyr show 
that significant cooling is going on at the cluster centre. The 
clusters are caught in the midst of ongoing feedback oscilla- 
tions, and because these oscillations are not in phase with 
each other, Tew(R) and S(R) at these small radii vary sig- 
nificantly across runs. 

Outside of 50 kpc, more subtle changes in temperature 
and entropy can be seen. In the left panels, in which the 
geometry of the bubbles is varied, the D5-8 run leads to 
slightly lower temperature and entropy profiles, the fiducial 
D5-10 stays close to the original profile, and the D5-12 run 
leads to a slight increase in temperature and entropy. The 
cluster as a whole is finding a slightly different equilibrium, 
which is weakly dependent on the geometry of energy input. 
On the other hand, in the right panels, in which the feedback 
efficiency is varied, all runs maintain Tew(R) and S(R) pro- 
files which are indistinguishable from each other. Again this 
similarity is indicative of self-regulation and occurs in spite 
of an order of magnitude variation in feedback efficiency. 

By t = 1.5 Gyr, the central oscillations have died down 
and the central entropies have been raised to values ^ 10 
keV cm 2 . Outside of this core region, the profiles have 
evolved slightly in the D5-8 and D5-12 cases, and are very 
similar to the initial state for all runs with rbubbie = 10 
kpc. As is observed in cool-core clusters, the AGN has man- 
aged to maintain the original entropy propfile in the face 
of widespread cooling. As we start already with a cluster 
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Figure 8. X-ray maps (projections of K{T)riHn e ) at 1 Gyr, projected perpendicular to the bubble axis. Top from left to right: D5-10, 
D5-8, D5-12. Bottom from left to right: D2-10, D10-10, D20-10. 



profile that has non-zero core entropy, we cannot say much 
about the origin of the e ntropy floor in non-cool-core clusters 
l|Cavagnolo et al.1 12009I ). However, these simulations, that 
cover a long physical time and many AGN outbursts, sug- 
gest that it is difficult to boost the core entropy to values 
~ 100 keV cm 2 with a string of relatively weak outbursts. 

Finally, in Fig. [12] we plot the metallicity in our 
simulated clusters. While previous simulations have mod- 
eled the effect of AGN activity on metals profiles {e.g. 
iTornatore et alJ 12004k iBriiggenl [20021 ). none of these have 
been self-regulating. As a result of the intermittent AGN ac- 
tivity, the metals are distributed over a large volume. Both 
metal density and metallicity are plotted in Fig. 1121 which 
shows that the AGN is effective at transporting metals. 

A secondary maximum develops in the metallicity, 
which moves outward at around 20 km/s. As a result, the 
profiles of the metallicity at times 1200, 1600 and 2000 Myr 
show dips near the centre. This might explain the abundance 
dips in the iron distribution that ha ve been observed in the 
Perseus and Centau rus clusters {e.g. iMatsushita et ai]|2007k 
ISanders et al 1 l2004k and for a discussion on metal abun- 
dance patterns in clusters also see Simionescu et al. 2008). 
In Perseus, the iron abundance dips from above 0.6 solar 
at around 40 kpc to below 0.5 solar near 15 kpc from the 
centre, and these numbers are similar to what we find here. 
However, unlike Perseus, the metal abundance in our sim- 



ulations increases again in the very centre of the cluster, 
although this would not occur if the metal p roduction had 
slowed down over time (jRenzini et al.l[l993h . Future simu- 
lations with a time-dependent metal production rate might 
shed further light on the origin of metallicity dips in some 
galaxy clusters. 

The profiles of the turbulent parameters, L and K, 
which we have not plotted here, look very similar across 
runs. Over most of the cluster, the turbulent velocities are 
^ 100 km s _1 , and the turbulent kinetic energy peaks at a 
radius w 20 kpc, at the outer edge of the bubbles where the 
turbulence is driven. The magnitude of turbulent velocities 
is in line with studies of resonant scattering in the ICM that 
find that the isotropic turbulent velocities on spatial scales 
smaller than 1 kpc are less than 100 km/s (Werner et al. 
2009). 

The effective turbulent diffusivities produced by the 
intermittent AGN ar e of the order of 50 kpc km/s 
l|Rebusco et al.ll2005h . iDavid fe Nulsenl [j2008r i have probed 
the Fe distribution in eight cool-core clusters and find that 
the Fe is significantly more extended than the stellar mass 
in the central galaxy. They find that diffusion coefficients of 
300-3000 kpc km/s are required to produce the observed Fe 
profiles. Based on these inferred diffusion coefficients, they 
conclude that heating by turbulent diffusion of entropy and 
turbulent dissipation can balance the radiative losses in their 
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Figure 9. X-ray maps (projections of A(T)njyn e ) at 1.5 Gyr, projected perpendicular to the bubble axis. Top from left to right: D5-10, 
D5-8, D5-12. Bottom from left to right: D2-10, D10-10, D20-10. 



sample. Our simulations confirm this. In real clusters, other 
factors such as mergers may contribute to the broadening 
of the metallicity distribution. However, the magnitude of 
metal transport in galaxy clusters can be explained by the 
action of AGN alone. 



4 DISCUSSION 

In this study, we have developed a complete and self- 
consistent three-dimensional model of AGN self-regulation 
in cool core clusters. In this picture cooling is balanced by 
periodic energy input from buoyant radio bubbles, whose 
turbulent evolution couples them to the inflowing cool gas, 
through mixing on scales that presently can only be cap- 
tured with subgrid modeling. When turbulence is property 
accounted for, our model manages to regulate itself for a 
range of geometries and feedback efficiencies. A phase of 
episodic heating is followed by a phase in which the AGN 
is quiet as the cool core reforms, which is followed in turn 
by another phase of episodic heating. Turbulence also helps 
to isotropise the injected energy over a turbulent timescale 
and remedies the problem o f channel formation reported by 
IVernaleo fc Reynolds! (|2006h . Our results reproduce a num- 
ber of key features seen in X-ray observations including: 

• temperature, density and entropy profiles that stay 



fairly monotonous over time, despite rapid cooling near the 
cluster centre; 

• X-rays profiles that are strongly peaked towards the 
core, and show clear evidence of large bubble cavities rising 
towards the cluster outskirts; and 

• metals that are well mixed into the outer regions of 
the cluster and exhibiting metallicity "dips" near the cluster 
centre. 

Despite these successes, however, there are a number 
of areas in which our model remains incomplete. In our 
approach, feedback is controlled by the total mass inflow 
through a sphere of r = 3 kpc, ignoring the details of 
the inner accretion flow, which will go through various 
phases as revealed by optical observations of filaments (e.g. 
iFabian et alJ I2008T ). Furthermore, our simulations do not 
include radiative heating, star formation, thermal conduc- 
tion, cosmic rays, or magnetic fields. However, it is unclear 
how important these processes are for AGN self-regulation. 
While thermal conduction may be important in establish- 
ing th e self-similar entropy profiles found outside cluster 
cores (jCavagnolo et al . 2009), its role in the cool cores of 
clusters is uncertain. Similarly, while cosmic rays are un- 
likely to supply much pressure support in galaxy clusters, 
cosmic-ray pressure gradients lead to convection, which in 
turn heats the ICM by advecting internal energy inward 
and allowing the cosmic rays to do pdV work on the ther- 
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Figure 10. Profiles of density (top), entropy (middle) and emission-weighted temperature (bottom) at 1.0 Gyr. In all panels the fat 
solid line shows the initial profile. Left panels: solid: D5-10; dotted: D5-12 and short dashed: D5-8. Right panels: long dashed: D10-10; 
dot-short dashed: D2-10 and dot-long dashed: D20-10. 



mal plasma (jChandran fc Ras cra 2007). Recently, there has 
also been renewed interest in the effect of magnetic fields 
on the dynamics of the ICM, which is subject to a num- 
ber of instabil i ties including the magneto thermal instability 
l|Balbusll2000l : IParrish fc Quataertl 120081) and th e heat- flux 
driven buoyancy instability ( Parrish et al.l [ioOct ) . These in- 
stabilities can lead to large-scale convective motions under 
a wide range of condit ions. However as d emonstrated by 
numerical simulations l|Parrish et al.l 120081 ). these motions 
tend to shut themselves off and thus do not play a major 
role in cool-core heating. Nonetheless, in future simulations, 
each of these effects should be explored. 

The evolution of AGN in clusters involves an enormous 
range of time-scales, from the few days characteristic of the 



dynamical time at the event horizon to the few gigayears 
required for the development of a cooling catastrophe in the 
cluster core. Here we advance a new explanation for the 
physics that determines the frequency of the intermittent 
bubbles that regulate cool-core clusters. Turbulence is cru- 
cial for mediating the heat input by AGN to the region that 
fuels the supermassive black hole, and thus the time scale for 
the process is determined by the turbulent turn-over time. 
This time scale, which is around a few 10 7 years in our sim- 
ulations, is given by the cluster scale radius divided by the 
central sound speed, and thus depends on overall cluster 
properties rather than the details of energy injection. 

A number of observational methods have been employed 
to quantify how long an average radio source spends in 
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Figure 11. Profiles of density (top), entropy (middle) and emission-weighted temperature (bottom) at 1.5 Gyr. In all panels the fat 
solid line shows the initial profile. Left panels: solid: D5-10; dotted: D5-12 and short dashed: D5-8. Right panels: long dashed: D10-10; 
dot-short dashed: D2-10 and dot-long dashed: D20-10. 



an active state and the time between active phases. One 
method inv okes spectral ageing and lobe expansion speed 
arguments I Alexander fc Leahvl Il987h . yielding times of a 
few 10 7 years. Energy injection rates required to quench 
cooling flows p oint to time-scales of the order of few 10 7 to 
10 s years (e.g. lOwen fc Eilekll 19981 ; iMcNamara et al.ll2005l ; 



iNulsen et alJl2005h . Observations of multiple cavities in the 
Perseus and Virgo clusters also suggest times « 10 8 years 
between subsequent outbursts of the radio-loud AGN. Low- 
frequency radio surveys in clusters might also reveal radio 
ghosts whose synchrotron ages will provide important clues 
about the histories of self-regulating radio- loud AGN. 
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Figure 12. Left: metal density of our fiducial run D5-10 at times t=400, 800, 1200, 1600 and 2000 Myr. The sold black line shows the 
expected metal density after 2 Gyr in the absence of any AGN activity. Right: The corresponding metallicity in units of solar metallicity 
(50 times metal mass fraction). Note that times 1200, 1600 and 2000 Myr show metallicity dips near the centre. 
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